program main

!Este programa resolve um problema anisotr�pico com fonte exponencial
! d/dx I(x) - A I(x) = F e^(-x/uo)
! A � gerado em Gera_A
! O vetor F � gerado em gera_Fonte (parte constante da fonte)


! Obs: Este programa est� adaptado para calcular as componentes de Fourier
! do problema sem simetria azimutal, assim quando m=0 o fluxo � dividido por 2
! No caso com simetria azimutal, fazer m=0, mas tirar esta divis�o.

!use portlib
!use msimsl

implicit none

include './rotinascomuns/include/ltsnparameters-lv173.f90'

double precision   :: z(R+1),Chl(R),ap(R),bp(R),cp(R)
double precision   :: part1,part2,cte,a_w(NM),a_c(NM),wl(NM)

! DGEEQU parameters
real(kind=8)                        ::Rr(R*no),Cc(R*no)
real(kind=8)                        ::ROWCND,COLCND,AMAX
integer                              ::INFO
! DGEEQU parameters

integer               :: k,kk,iwl,pos,ptoss !,div,ptos
real(kind=8), dimension(div)      :: abcissa_pos
real(kind=8), dimension(div)      :: abcissa_neg
real(kind=8), dimension(div)      :: ordenada_pos
real(kind=8), dimension(div)      :: ordenada_neg
real(kind=8), dimension(ptos)      :: saida_xpos
real(kind=8), dimension(ptos)      :: saida_xneg
real(kind=8), dimension(ptos)      :: saida_ypos
real(kind=8), dimension(ptos)      :: saida_yneg
real(kind=8), dimension(2*ptos,R+1)   :: spr_valo_flu

 character(len=18)   :: anidir
 character(len=7)   :: profile,fnametemp
 character(len=18)       :: radfile
 character(len=2)        :: iwlchar
 character(len=200)   :: fileparam,filedata
 character(len=200)   :: fileprofile,filewl

! character(len=17)   :: anidir
! character(len=50)   :: fileparam,filedata,filerad

include './rotinascomuns/include/pimparameters.f90'

!anidir = '../dados/entrada/'
!fileparam='../dados/entrada/parametros_hidro_10r_lv0_pnf'
!filedata = '../dados/saida/50/data_omega_pnf.dat'
!filerad = '../dados/saida/50/radiance_pnf.dat'


anidir = './dados/clorofila/'

 call getarg(1,profile) 
 call getarg(2,radfile)   
 call getarg(3,iwlchar)   

fnametemp = 'temp.dat'
write (fnametemp,*) iwlchar
read (fnametemp,*) iwl

fileparam = anidir//'perfis/'//profile//'/parametros_hidro_'//profile//'.dat'
fileprofile = anidir//'perfis/'//profile//'/'//profile//'.dat'
filewl = anidir//'eq327_m11.dat'
!filedata = anidir//'radiancias/'//profile//'/radiancias.dat'
filedata = anidir//'radiancias/'//profile//'/'//radfile

!dados do problema

OPEN (21, FILE = fileparam)
READ (21, *) j
READ (21, *) cc0(1), ccL(1), valor_tau
READ (21, *) uo
READ (21, *) ptoss

READ (21, *) (saida_xpos(i), i = 1, ptos)
DO i = 1, ptos
   saida_xneg(i) = - saida_xpos(ptos+1-i)
ENDDO

OPEN (21, FILE = fileprofile)
DO i = 1, max(R+1,2)
   READ (21,*) z(i),Chl(i)
      WRITE(*,"(2E13.5)") z(i),Chl(i)
ENDDO
CLOSE(21)
         
OPEN (21, FILE = filewl)
DO i = 1, NM
   READ (21, *) wl(i),a_w(i),a_c(i)
ENDDO
CLOSE(21)


ipath=1

pi=3.141593d+00   
if (j==1) then
   phi=0.0d+00
elseif(j==2) then
   phi=pi/2.0d+00
else
   phi=pi
endif
!uo=0.5d+00                        ! expoente da fonte e^(-x/uo)
!open (unit=20, file='saida_hidro')         ! saida da respota sai.dat
!write(20,*) 'N=',no,'Grau de anisotropia, L =',lv+1
!write(20,*) 'm=',m
!Alocando Matrizes e Vetores 
!Definindo as CContorno nos vetores de ordem N/2 cc0 e ccL

cc0 = 0.0d+00
ccL = 0.0d+00



if(lv==0) then
   open (unit=8,file=anidir//'legendre/hl_0.dat',status='old')
elseif(lv==8) then
   open (unit=8,file=anidir//'legendre/hl_8.dat',status='old')      !arquivo dos Bl da anisotropia, B0=1
elseif(lv==82) then
   open (unit=8,file=anidir//'legendre/hl_82.dat',status='old')
elseif(lv==173) then
   open (unit=8,file=anidir//'legendre/hl_173.dat',status='old')
elseif(lv==299) then
   open (unit=8,file=anidir//'legendre/hl_299.dat',status='old')
endif
do i=1,lv+1
   read(8,*) hl(i)
enddo
close(8)

call galeg(no,-1.d+00,1.d+00,u,w)            ! calcula pesos e raizes da quadratura

u(no+1)=uo
!u(no+2)=0.0
kk=iwl
   part2 = 1 + 0.2*dexp(-0.014*(wl(kk) - 440.0))
   cte = (550.0/wl(kk)) * 0.30
   DO i = 1, R
      part1 = a_w(kk) + 0.06*a_c(kk)*(Chl(i)**0.65)
      ap(i) = part1*part2
      bp(i) = cte*(Chl(i)**0.62)
      cp(i) = ap(i) + bp(i)
      omega(i) = bp(i) / cp(i)
      L(i) = cp(i)*(z(i+1)-z(i))
      WRITE(*,"(5E13.5)") ap(i),bp(i),cp(i),omega(i),L(i)
   ENDDO

   LL(1)=0.0d+00
   IF (R .GT. 1) THEN
      do i=2,R+1
         LL(i)=LL(i-1)+L(i-1)
      enddo
   ELSE
     LL(2) = L(1)
   ENDIF


soma0=0
somaL=0

DO MM=1,LV+1
   M=MM-1
   MatrizSist=0.0
   call pnass(m,u,lv,no,pol)                  ! calcula os polinomios associados de Legendre
   call geraA(m,no,A,hl,omega(1),lv,pol,u,w)         ! gera a matriz A de sI-A
   call gera_Fonte(m,no,hl,omega(1),lv,pol,u,fonte)   ! Adaptar a cada caso, gera o vetor cte da fonte
   call calc_dados(no,A,eigval,eigvec,inv_eigvec)   ! Autovalores, Autovetores, Inversa dos Autovetores
   call calc_b(no,eigval,eigvec,0.0d+00,L(1),b0)      ! B(0)=X Exp[x* lambda]
   call convolucao(no,eigval,eigvec,inv_eigvec,fonte,uo,0.0d+00,L(1),Hh0)
   do i=1,no/2    
      do j=1,no
         MatrizSist(i,j)=b0(i,j)
      enddo
   enddo
   do i=1,no/2
      VetorSist(i)=cc0(i)-Hh0(i)
   enddo
   do i=1,no
      do j=1,no
         Grande_B0(1,i,j)=b0(i,j)
      enddo
   enddo
   do i=1,no
      Grande_H0(1,i)=Hh0(i)
   enddo
   
   do jj=2,R
      !omega(jj)=0.1*(21./(dfloat(jj)+20.))**4
      call calc_b(no,eigval,eigvec,L(jj-1),L(jj-1),bL)            ! BL(L) e BNeg(L)
      call convolucao(no,eigval,eigvec,inv_eigvec,fonte,uo,LL(jj),LL(jj),HhL)
      do ij=1,no
         do ji=1,no
            MatrizSist(no/2+(jj-2)*no+ij,(jj-2)*no+ji)=bL(ij,ji)
         enddo
      enddo
      do i=1,no
         do j=1,no
            Grande_BL(jj-1,i,j)=bL(i,j)
         enddo
      enddo
      do i=1,no
         Grande_HL(jj-1,i)=HhL(i)
      enddo

      call geraA(m,no,A,hl,omega(jj),lv,pol,u,w)         ! gera a matriz A de sI-A
      call gera_Fonte(m,no,hl,omega(jj),lv,pol,u,fonte)   ! Adaptar a cada caso, gera o vetor cte da fonte
      call calc_dados(no,A,eigval,eigvec,inv_eigvec)   ! Autovalores, Autovetores, Inversa dos Autovetores
      call calc_b(no,eigval,eigvec,0.0d+00,L(jj),b0)      ! B(0)=X Exp[x* lambda]
      call convolucao(no,eigval,eigvec,inv_eigvec,fonte,uo,LL(jj),LL(jj+1),Hh0)
      do ij=1,no
         do ji=1,no
            MatrizSist(no/2+(jj-2)*no+ij,(jj-1)*no+ji)=-b0(ij,ji)
         enddo
      enddo
      do ii=1,no
         VetorSist(no/2+(jj-2)*no+ii)=Hh0(ii)-HhL(ii)
      enddo
      do i=1,no
         do j=1,no
            Grande_B0(jj,i,j)=b0(i,j)
         enddo
      enddo
      do i=1,no
         Grande_H0(jj,i)=Hh0(i)
      enddo
      !print*, jj
   enddo
   call calc_b(no,eigval,eigvec,L(R),L(R),bL)            ! BL(L) e BNeg(L)
   call convolucao(no,eigval,eigvec,inv_eigvec,fonte,uo,LL(R+1),LL(R+1),HhL)
   do ij=1,no/2
      do ji=1,no
         MatrizSist(R*no-no/2+ij,R*no-no+ji)=bL(ij+no/2,ji)
      enddo
   enddo
   do ii=1,no/2
      VetorSist(R*no-no/2+ii)=ccL(ii+no/2)-HhL(ii+no/2)
   enddo
   do i=1,no
      do j=1,no
         Grande_BL(R,i,j)=bL(i,j)
      enddo
   enddo
   do i=1,no
      Grande_HL(R,i)=HhL(i)
   enddo


   CALL DGEEQU( R*no, R*no, MatrizSist, R*no, Rr, Cc, ROWCND, COLCND, AMAX,INFO)
   !write(*,*) 'INFO DGEEQU: ',INFO
   do ij=1,R*no
      do ji=1,R*no
         MatrizSist(ij,ji) = Rr(ij)*MatrizSist(ij,ji)*Cc(ji)
      enddo
   enddo

    do ij=1,R*no
      VetorSist(ij) = Rr(ij)*VetorSist(ij)
   enddo

! Compute preconditioners
      IF (PRET.EQ.1) THEN
          DO I = 1,N
              Q1(I) = 1.0D0/MatrizSist(I,I)
          ENDDO
      ELSE IF (PRET.EQ.2) THEN
          DO I = 1,N
              Q2(I) = 1.0D0/MatrizSist(I,I)
          ENDDO
      ELSE IF (PRET.EQ.3) THEN
          DO I = 1,N
              Q1(I) = 1.0D0/SQRT(MatrizSist(I,I))
              Q2(I) = Q1(I)
          ENDDO
      END IF
! Set polynomial preconditioner coefficients
      IF (POLYT.EQ.1) THEN
          DO I = 1,MMM + 1
              GAMMA(I) = 1.0D0
          ENDDO
      ELSE IF (POLYT.EQ.2) THEN
          CALL DWLSCOEF(MMM,GAMMA)
      ELSE IF (POLYT.EQ.3) THEN
          CALL DULSCOEF(MMM,GAMMA)
      END IF
    
    !Rotinas PIM

    SaiSist = 0.0D0
    CALL PIMDSETPAR(IPAR,DPAR,LDA,N,N,N,C,-1,-1,PRET,STOPT,MAXIT,TOL)
   CALL PIMDCGNR(SaiSist,VetorSist,WRK,IPAR,DPAR,MATVEC,TMATVEC,POLYL,POLYR,PDSUM,PDNRM2,PROGRESS)

    do ij=1,R*no
      SaiSist(ij) = SaiSist(ij)*Cc(ij)
   enddo

   !write(*,*) 'm, ipar(12), ipar(13) =',m,ipar(12),ipar(13)

    do ij=1,no
       do ji=1,no
         b0(ij,ji)=Grande_b0(1,ij,ji)
      enddo
      Hh0(ij)=Grande_H0(1,ij)
    enddo
   SaiSist_Aux=SaiSist(1:no)
   call Fluxo_angular(no,SaiSist_Aux,b0,Hh0,Fluxo_em_0)
      do ii=1,no
      fluxo_em_0(ii)=fluxo_em_0(ii)*dcos(m*phi)
    enddo
    soma0(1,:)=soma0(1,:) + fluxo_em_0

    do i=1,R-1
      do ij=1,no
         do ji=1,no
            bL(ij,ji)=Grande_bL(i,ij,ji)
            b0(ij,ji)=Grande_b0(i+1,ij,ji)
         enddo
         HhL(ij)=Grande_HL(i,ij)
         Hh0(ij)=Grande_H0(i+1,ij)
      enddo
      SaiSist_Aux=SaiSist((i-1)*no+1:i*no)
       call Fluxo_angular(no,SaiSist_Aux,bL,HhL,Fluxo_em_L)
      SaiSist_Aux=SaiSist(i*no+1:(i+1)*no)
       call Fluxo_angular(no,SaiSist_Aux,b0,Hh0,Fluxo_em_0)
       do ii=1,no
         fluxo_em_L(ii)=fluxo_em_L(ii)*dcos(m*phi)
         fluxo_em_0(ii)=fluxo_em_0(ii)*dcos(m*phi)
       enddo
       somaL(i,:) = somaL(i,:) + fluxo_em_L
      soma0(i+1,:) = soma0(i+1,:) + fluxo_em_0
   enddo

    do ij=1,no
       do ji=1,no
         bL(ij,ji)=Grande_bL(R,ij,ji)
      enddo
      HhL(ij)=Grande_HL(R,ij)
    enddo
   fluxo_em_L=0
   SaiSist_Aux=SaiSist(R*no-no+1:R*no)
   call Fluxo_angular(no,SaiSist_Aux,bL,HhL,Fluxo_em_L)
   !fluxo_em_L(no/2+1:no) = Matmul(bL(no/2+1:no,:),Sai_Gauss_Aux) + HhL(no/2+1:no)
   !fluxo_em_L(1:no/2) = Matmul(bL(1:no/2,:),Sai_Gauss_Aux) + HhL(1:no/2)
      do ii=1,no
      fluxo_em_L(ii)=fluxo_em_L(ii)*dcos(m*phi)
    enddo
    somaL(R,:) = somaL(R,:) + fluxo_em_L
   !write(*,*) 'm, regiao, soma : ',m ,R, sum(fluxo_em_L(no/2+1:no))
enddo


!----------------------------------------------------------------
! Aqui eu chamo os mus continuos, via interpolacao 433
! Detalhe importante eh que eu tenho que dividir tudo em
! arrays auxiliares para mandar para dentro da rotina. Na
! volta eu junto tudo de novo
!----------------------------------------------------------------


do j=1,div
   abcissa_pos(j)=u(div-j+1)
   abcissa_neg(j)=u(no-j+1)
end do

do k=1,R
   do j=1,div
      pos=div-j+1
      ordenada_pos(j)=soma0(k,pos)
      ordenada_neg(j)=soma0(k,no-j+1)
   end do
   call intrpl(2,div,abcissa_pos,ordenada_pos,ptos,saida_xpos,saida_ypos)
   call intrpl(2,div,abcissa_neg,ordenada_neg,ptos,saida_xneg,saida_yneg)
   do j=1,ptos
      spr_valo_flu(j,k)=saida_yneg(j)
      spr_valo_flu(j+ptos,k)=saida_ypos(j)
   end do
end do

   do j=1,div
      pos=div-j+1
      ordenada_pos(j)=somaL(R,pos)
      ordenada_neg(j)=somaL(R,no-j+1)
   end do
   call intrpl(2,div,abcissa_pos,ordenada_pos,ptos,saida_xpos,saida_ypos)
   call intrpl(2,div,abcissa_neg,ordenada_neg,ptos,saida_xneg,saida_yneg)
   do j=1,ptos
      spr_valo_flu(j,R+1)=saida_yneg(j)
      spr_valo_flu(j+ptos,R+1)=saida_ypos(j)
   end do


  OPEN(22,FILE=filedata,STATUS="UNKNOWN", &
               ACCESS="SEQUENTIAL",ACTION="WRITE")
  DO i = 1, ptos
     WRITE (22, "(1P40E11.3)") (spr_valo_flu(i,j), j=1,R+1)
  ENDDO
  DO i = ptos+1,2*ptos
     WRITE (22, "(1P40E11.3)") (spr_valo_flu(i,j), j=1,R+1)
  ENDDO
  CLOSE(22)


end
